/*---------- Ｍａｉｎ ---------*/

main() {

double t, t1, t3, tz=0 , dt=TO/NM ,irw=ir_w ;

double p, k , w0=w*ir_w, h_w=HH*ir_h ;

double z, z0=0, zmax=0 ;

double v, ,0=0.1 ;

double C, D , D1, D2 ;

double k1, k2, k3, k4 ;

double L1, L2, L3, L4 ;

double n=0 flg=0 ;

double prd1, prd2 ;

double h0=h0, dc_=dc ;

/**/

k = kn () ;

D1= h_w * 0.75 * w0 * cosh(k*(h+HWL-dc)) / cosh(k*(h+HWL)) ;

D2= h_w * 0.5 * w0 * cosh(k*(h+HWL-dc)) / cosh(k*(h+HWL)) ;

/**/

if (h0_>0.75*h_w) t1=TO

else t1=asin(h0_/0.75/h_w)/w0 ;

if (dc_>0.5*h_w) t3=TO*irw ;

else t3=asin(dc_/0.5/h_w) /w0 + TO/(2*irw) ;

v = v0 ; z = z0 ;

fprintf(stdprn,"%s\n",NAME) ;

fprintf(stdprn,"波 高＝%4.1f 波 数＝ %6.4\n",HH,k ) ;

fprintf(stdprn,"周 期＝%4.1f 角速度＝%6.3f\n",TO,w ) ;

fprintf(stdprn,"水 深＝%4.1f 圧力P'＝%6.3f\n",h ,D1) ;

fprintf(stdprn,"天端高＝%4.1f H.W.L ＝%4.1f\n",h0_+HWL,HWL) ;

fprintf(stdprn,"没水部＝%4.1f\n",dc_-HWL) ;

fprintf(stdprn,\n 時刻 圧力 位置 速度 加速度\n");

/**/

for (t=t1 ; t<=t3 ; t+=dt) {

if ( fmod(floor(2*irw*t1/TO) ,2) == 0) D=D1 ;

else D=D2 ;

p = D*sin(w0*t) - h0_;

if (v>=0) C=C1 ; /* 上昇時係数 */

else C=C2 ; /* 下降時係数 */

k1 = alfa(t , v , C, D)*dt ; L1 = v *dt ;

k2 = alfa(t+dt/2, v+k2/2 , C, D)*dt ; L2 = (v+k1/2)*dt ;

k3 = alfa(t+dt/2, v+k2/2 , C, D)*dt ; L3 = (v+k2/2)*dt ;

k4 = alfa(t+dt , v+k3 , C, D)*dt ; L4 = (v+k3) *dt ;

/**/

v += (k1+2*k2+2*k3+k4)/6 ; z += (L1+2*L2+2*L3+L4)/6 ;

if (fabs(v)<0.01) v=fabs(v)/v*0.01 ;

/**/

if (z>zmax) zmax = z ;

if (z<0 && flg==0) { tz=t; flg=1 ;}

/**/

if (fmod(n, (NM/40))==0)

fprintf(stdprn,"%8.3f %8.3f %8.3f %8.3f %8.3f\n",t,p,z,v,k1/dt) ;

n += 1 ;

}

fprintf(stdprn,"%8.3f %8.3f %8.3f %8.3f %8.3f\n",t,p,z,v,k1/dt) ;

fprintf(stdprn,"\n Zmax=%6.3f t(z=0)=%5.3f\n",zmax,tz) ;

/**/

if (z>0) {

prd1=1.61*sqrt(z) ; prd2=To/ir_w+t1-t3 ;

fprintf(stdprn,"自由落下時間=%6.1f 再侵入迄の時間%6.1f\n",prd1,prd2) ;

}

}

/*---------- Program end ----------*/

前ページ　　　目次へ　　　次ページ